In [1]:
import scipy.integrate as spi
import numpy as np
import pylab as pl
from matplotlib.pyplot import figure
import pandas as pd
In [2]:
import requests
import urllib.request
import time
from bs4 import BeautifulSoup
In [3]:
#Link to install platypus
# https://github.com/Project-Platypus/Platypus?fbclid=IwAR3VssR7uHGCRiROvivIbeUwUboBWLaO3qSMvVwlH7Rja2nh4rBrLheXyeI
#Genetic Algorithm Library
from platypus.algorithms import NSGAII, Problem
#from platypus.algorithms import NSGAII, Problem, Real
from platypus.config import Real
In [ ]:
# Update the dataset (save the updated data)
import urllib.request
print('Beginning file download with urllib2...')
url = 'https://raw.githubusercontent.com/wcota/covid19br/master/cases-brazil-states.csv'
urllib.request.urlretrieve(url, '/media/juliane_oliveira/My Passport/dropbox_15_11_2019/Artigos_preprints/COVID19/data_brasil/estadosnet.csv')
In [4]:
# run code to organize data
import datetime as dt
def ler_banco(arq,var):
banco = pd.read_csv(arq)
banco =banco[banco[var].notnull()]
if var=='cod_city':
banco[var] = pd.to_numeric(banco[var],downcast='integer')
nome_local =list(banco[var].unique())
for i in banco.index:
banco.date[i] = dt.datetime.strptime(banco.date[i], '%Y-%m-%d').date()
nome_local =list(banco[var].unique())
local = []
for est in nome_local:
aux = banco[banco[var]==est].sort_values('date')
data_ini = aux.date.iloc[0]
data_fim = aux.date.iloc[-1]
dias = (data_fim-data_ini).days + 1
d = [(data_ini + dt.timedelta(di)) for di in range(dias)]
if var=='cod_city':
cod_city = [est for di in range(dias)]
estado = [aux.state.iloc[0] for di in range(dias)]
uf = [aux.UF.iloc[0] for di in range(dias)]
city = [aux.city.iloc[0] for di in range(dias)]
df = pd.DataFrame({'date':d,'state':estado,'UF':uf,'city':city,var:cod_city})
df.cod_city =pd.to_numeric(df.cod_city,downcast='integer')
else:
estado = [est for di in range(dias)]
df = pd.DataFrame({'date':d,var:estado})
casos = []
caso = 0
i_aux = 0
for i in range(dias):
if (d[i]-aux.date.iloc[i_aux]).days==0:
caso = aux.totalCases.iloc[i_aux]
casos.append(caso)
i_aux=i_aux+1
else:
casos.append(caso)
df['totalCases'] = casos
local.append(df)
return nome_local, local
In [5]:
# read the updated data you downloaded in cell 4
df= ler_banco('/media/juliane_oliveira/My Passport/dropbox_15_11_2019/Artigos_preprints/COVID19/data_brasil/estadosnet.csv','state')
In [6]:
pop = pd.read_excel('/media/juliane_oliveira/My Passport/dropbox_15_11_2019/Artigos_preprints/COVID19/data_brasil/popUF.xlsx')
In [7]:
#n = 27
for n in range(len(pop)):
UF = df[0][n] #df.state.unique()[n]
dfUF = np.array(df[1][n]['totalCases'])
date = df[1][n]['date'][0]
print(date)
popUF = pop[pop.UF == UF]['pop'].unique()
dfUF_norm = dfUF/popUF
I0 = dfUF_norm[1]
S0= 1 -I0
#print('condições iniciais para ' + UF)
data = dfUF_norm
pop2 = popUF
def diff_eqs(INP,t,beta,gamma):
'''The main set of equations'''
Y=np.zeros((3))
V = INP
Y[0] = - beta * V[0] * V[1]
Y[1] = beta * V[0] * V[1] - gamma * V[1]
Y[2] = gamma * V[1]
return Y # For odeint
#Optimisation Function which minimises the MSE between the SIR model and the data
def fitness_function(x):
beta = x[0]
gamma = x[1]
result = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
mean_squared_error = ((np.array(data)-result[:,1])**2).mean()
return [mean_squared_error]
#GA Parameters
number_of_generations = 10000
ga_population_size = 100
number_of_objective_targets = 1 # The MSE
number_of_constraints = 0
number_of_input_variables = 2 # beta and gamma
problem = Problem(number_of_input_variables,number_of_objective_targets,number_of_constraints)
problem.types[0] = Real(0, 1) #beta initial Range
print('Resultados gerais para o estado ' + UF)
print(' Gamma ' + ' Beta ' + ' R0 ')
for j in range(3):
if (j == 0):
problem.types[1] = Real(1/5-0.00000001, 1/22+-0.00000001) #gamma initial Range
#problem.types[1] = Real(1/14-0.00000001, 1/20) #
#problem.types[1] = Real(1/20-0.00000001, 1/22+-0.00000001)
problem.function = fitness_function
algorithm = NSGAII(problem, population_size = ga_population_size)
#The model needs to output in the same shape as the data to be fitted
TS=1
ND=len(data)-1
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
INPUT = (S0, I0, 0.0) # R(0) = 0
#Running the GA
algorithm.run(number_of_generations)
#Getting the feasible solutions only just to be safe
#The GA may generate bad solutions if number_of_generations is not enough or there are too many contraints
feasible_solutions = [s for s in algorithm.result if s.feasible]
beta = feasible_solutions[0].variables[0]
#print('Beta Optimsed = {}'.format(round(beta,2)))
gamma= feasible_solutions[0].variables[1]
#print('Gamma Optimsed = {}'.format(round(gamma,2)))
R0 = beta/gamma
#print('R0 = {}'.format(round(R0,2)))
print(' 1/5 ' + ' {} '.format(round(beta,2)) + ' {}'.format(round(R0,2)))
#Running the model and plotting with the optimised values
ND=2*(len(data)-1)
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
pl.plot(data*100, c = 'red', label=UF,linewidth=5)
pl.plot(RES[:,1]*100, label='SIR Model')
pl.xlabel('Dias',fontsize=20)
pl.ylabel('Casos Confirmados (%)',fontsize=20)
pl.xticks(range(1,t_end+1,8),fontsize=10)
pl.grid(alpha = 0.5,which='both')
pl.legend(fontsize=15)
pl.show()
pl.semilogy(RES[:,1]*pop2, '-r', label='Infectados')
pl.semilogy(data*pop2, 'ob')
#pl.loglog(time2,func(time2)/14873064, '-g')
pl.xlim(1,30)
pl.xlabel('Dias')
pl.ylabel('Infectados')
pl.show()
# Save projection
dataset = pd.DataFrame({'S': RES[:, 0]*pop2, 'I': RES[:, 1]*pop2, 'R':RES[:, 2]*pop2})
dataset.to_csv('predito26_03_{}_gamma1_14.csv'.format(UF))
#Now, running the model for a longer period
TS=0.1
ND=100
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
#Ploting of the fitted curve only
from matplotlib.pyplot import figure
figure(figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
pl.subplot(311)
pl.plot(t_range,RES[:,0], '-g', label='SuscetÃveis')
pl.plot(t_range,RES[:,2], '-k', label='Recuperados')
pl.plot(t_range,RES[:,1], '-r', label='Infectado')
pl.plot(data, 'ob')
pl.ylim(0,1)
pl.legend(loc=0)
pl.title('Dinâmica do CoviD19 {}: R0 = {}, beta = {}, gamma = {}. '.format(UF,round(R0,2),round(beta,3),round(gamma,3)))
pl.ylabel('Proporção da população')
pl.xlabel('Tempo (dias)')
pl.savefig('plot1_25_03_{}_gamma1_5.pdf'.format(UF))
pl.show()
elif (j == 1):
#problem.types[1] = Real(1/5-0.00000001, 1/22+-0.00000001) #gamma initial Range
problem.types[1] = Real(1/14-0.00000001, 1/20) #
#problem.types[1] = Real(1/20-0.00000001, 1/22+-0.00000001)
problem.function = fitness_function
algorithm = NSGAII(problem, population_size = ga_population_size)
#The model needs to output in the same shape as the data to be fitted
TS=1
ND=len(data)-1
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
INPUT = (S0, I0, 0.0) # R(0) = 0
#Running the GA
algorithm.run(number_of_generations)
#Getting the feasible solutions only just to be safe
#The GA may generate bad solutions if number_of_generations is not enough or there are too many contraints
feasible_solutions = [s for s in algorithm.result if s.feasible]
beta = feasible_solutions[0].variables[0]
#print('Beta Optimsed = {}'.format(round(beta,2)))
gamma= feasible_solutions[0].variables[1]
#print('Gamma Optimsed = {}'.format(round(gamma,2)))
R0 = beta/gamma
#print('R0 = {}'.format(round(R0,2)))
print(' 1/14 ' + ' {} '.format(round(beta,2)) + ' {}'.format(round(R0,2)))
#Running the model and plotting with the optimised values
ND=2*(len(data)-1)
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
pl.plot(data*100, c = 'red', label=UF,linewidth=5)
pl.plot(RES[:,1]*100, label='SIR Model')
pl.xlabel('Dias',fontsize=20)
pl.ylabel('Casos Confirmados (%)',fontsize=20)
pl.xticks(range(1,t_end+1,8),fontsize=10)
pl.grid(alpha = 0.5,which='both')
pl.legend(fontsize=15)
pl.show()
pl.semilogy(RES[:,1]*pop2, '-r', label='Infectados')
pl.semilogy(data*pop2, 'ob')
#pl.loglog(time2,func(time2)/14873064, '-g')
pl.xlim(1,30)
pl.xlabel('Dias')
pl.ylabel('Infectados')
pl.show()
# Save projection
dataset = pd.DataFrame({'S': RES[:, 0]*pop2, 'I': RES[:, 1]*pop2, 'R':RES[:, 2]*pop2})
dataset.to_csv('predito26_03_{}_gamma1_14.csv'.format(UF))
#Now, running the model for a longer period
TS=0.1
ND=100
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
#Ploting of the fitted curve only
from matplotlib.pyplot import figure
figure(figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
pl.subplot(311)
pl.plot(t_range,RES[:,0], '-g', label='SuscetÃveis')
pl.plot(t_range,RES[:,2], '-k', label='Recuperados')
pl.plot(t_range,RES[:,1], '-r', label='Infectado')
pl.plot(data, 'ob')
pl.ylim(0,1)
pl.legend(loc=0)
pl.title('Dinâmica do CoviD19 {}: R0 = {}, beta = {}, gamma = {}. '.format(UF,round(R0,2),round(beta,3),round(gamma,3)))
pl.ylabel('Proporção da população')
pl.xlabel('Tempo (dias)')
pl.savefig('plot1_25_03_{}_gamma1_14.pdf'.format(UF))
pl.show()
else:
#problem.types[1] = Real(1/5-0.00000001, 1/22+-0.00000001) #gamma initial Range
#problem.types[1] = Real(1/14-0.00000001, 1/20) #
problem.types[1] = Real(1/20-0.00000001, 1/22+-0.00000001)
problem.function = fitness_function
algorithm = NSGAII(problem, population_size = ga_population_size)
#The model needs to output in the same shape as the data to be fitted
TS=1
ND=len(data)-1
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
INPUT = (S0, I0, 0.0) # R(0) = 0
#Running the GA
algorithm.run(number_of_generations)
#Getting the feasible solutions only just to be safe
#The GA may generate bad solutions if number_of_generations is not enough or there are too many contraints
feasible_solutions = [s for s in algorithm.result if s.feasible]
beta = feasible_solutions[0].variables[0]
#print('Beta Optimsed = {}'.format(round(beta,2)))
gamma= feasible_solutions[0].variables[1]
#print('Gamma Optimsed = {}'.format(round(gamma,2)))
R0 = beta/gamma
#print('R0 = {}'.format(round(R0,2)))
print(' 1/21 ' + ' {} '.format(round(beta,2)) + ' {}'.format(round(R0,2)))
#Running the model and plotting with the optimised values
ND=2*(len(data)-1)
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
pl.plot(data*100, c = 'red', label=UF,linewidth=5)
pl.plot(RES[:,1]*100, label='SIR Model')
pl.xlabel('Dias',fontsize=20)
pl.ylabel('Casos Confirmados (%)',fontsize=20)
pl.xticks(range(1,t_end+1,8),fontsize=10)
pl.grid(alpha = 0.5,which='both')
pl.legend(fontsize=15)
pl.show()
pl.semilogy(RES[:,1]*pop2, '-r', label='Infectados')
pl.semilogy(data*pop2, 'ob')
#pl.loglog(time2,func(time2)/14873064, '-g')
pl.xlim(1,30)
pl.xlabel('Dias')
pl.ylabel('Infectados')
pl.show()
# Save projection
dataset = pd.DataFrame({'S': RES[:, 0]*pop2, 'I': RES[:, 1]*pop2, 'R':RES[:, 2]*pop2})
dataset.to_csv('predito26_03_{}_gamma1_21.csv'.format(UF))
#Now, running the model for a longer period
TS=0.1
ND=100
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range,args=(beta,gamma))
#Ploting of the fitted curve only
from matplotlib.pyplot import figure
figure(figsize=(10, 10), dpi=80, facecolor='w', edgecolor='k')
pl.subplot(311)
pl.plot(t_range,RES[:,0], '-g', label='SuscetÃveis')
pl.plot(t_range,RES[:,2], '-k', label='Recuperados')
pl.plot(t_range,RES[:,1], '-r', label='Infectado')
pl.plot(data, 'ob')
pl.ylim(0,1)
pl.legend(loc=0)
pl.title('Dinâmica do CoviD19 {}: R0 = {}, beta = {}, gamma = {}. '.format(UF,round(R0,2),round(beta,3),round(gamma,3)))
pl.ylabel('Proporção da população')
pl.xlabel('Tempo (dias)')
pl.savefig('plot1_25_03_{}_gamma1_21.pdf'.format(UF))
pl.show()
print('\r\n')
FIM